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In order to discover generic effects of heterogeneous communication delays on the dynamics of large 
systems of coupled oscillators, this paper studies a modification of the Kuramoto model incorporating 
a distribution of interaction delays. By focusing attention on the reduced dynamics on an invariant 
manifold of the original system, we derive governing equations for the system which we use to study 
stability of the incoherent states and the dynamical transitional behavior from stable incoherent 
states to stable coherent states. We find that spread in the distribution function of delays can 
greatly alter the system dynamics. 



PACS numbers: 05.45.Xt,05.45.-a,89.75.-k 



Large systems of coupled oscillators occur in many sit- 
uations in modern science and engineering [l| . Noted ex- 
amples include synchronous flashing of fireflies Q , pedes- 
trian induced oscillations of the Millennium Bridge 
cardiac pace-maker cells alpha rhythms in the brain 
Q, glycolytic oscillations in yeast populations Q, cellu- 
lar clocks governing circadian rhythm in mammals 0], 
oscillatory chemical reactions [§[, etc. An incompletely 
understood aspect of such systems is that signal prop- 
agation may take non-negligible time, and that systems 
often have a finite reaction time to inputs that they re- 
ceive. Time delays are thus both natural and inevitable 
in many of these systems. In order to elucidate phenom- 
ena induced by time delay in large coupled oscillator sys- 
tems, Refs.[9l.lirj| and carried out studies of globally 
coupled oscillators of the Kuramoto type [ijj] in the pres- 
ence of time delay. These previous works all treated the 
special case in which all time delays between interacting 
oscillators were identical, and, in that context, they un- 
covered many interesting behaviors revealing that time 
delay can profoundly affect the dynamics of coupled os- 
cillator systems. However, in most situations where de- 
lays are an important consideration, the delays are not 
all identical. The aim of this paper is to study the more 
realistic case where there is a distribution of time de- 
lays along the links connecting the oscillators. We shall 
see that previous striking features obtained in the case 
of uniform time delay are evidently strongly dependent 
on coherent communication between oscillators, and, as 
a consequence, are substantially changed by the incor- 
poration of even modest spread in the time delays. For 
example, comparing results for typical cases with uni- 
form delay and with a 30% spread in delay, we will show 
that this delay spread (a) can completely eliminate the 
resonant structure in the average delay time dependence 
of the critical coupling k c for the onset of coherence, (b) 
can introduce hysteresis into the system behavior, and 
(c) can substantially decrease the number of attractors 
that simultaneously exist in a given situation. 

We consider a network of oscillators with all-to-all cou- 
pling according to the classical Kuramoto scheme, but 
incorporating link-dependent interaction time delays Tij 



for coupling between any two oscillators i and j, 
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where 9i(t) is the phase of oscillator i, oji is the natu- 
ral frequency of oscillator i, k characterizes the coupling 
strength between oscillators, N is the total number of 
oscillators, tu = 0, and i = 1,2- --N. Following Ku- 
ramoto, we note that the effect of all the oscillators in 
the network on oscillator i may be expressed in terms of 
an "order parameter" r^, 
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To facilitate the analysis, we consider the following two 
simplifying assumptions. First, we consider the contin- 
uum limit N — * 00 appropriate to the study of large 
systems, N 3> 1. Second, we assume the collection of 
all delays is characterized by a distribution h{r) such 
that the fraction of links with delays between r and r+dr 
is h(r)dr. We, furthermore, assume that, for randomly 
chosen links, r is uncorrelated with the oscillator fre- 
quencies lu at cither end of the link. These assumptions 
enable a description of the system dynamics in terms of 
a single oscillator distribution function f(6,u>,t), which 
evolves in response to a mean field r(t) according to the 
following oscillator continuity equation, 
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In this case, the mean field r(t) is given by 



r(t) = / ^ - r)h(r)dr, 
Jo 



/ = 0. (4) 
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where Eq. ([6]) gives the input that nodes would receive in 
the absence of delay, and Eq.® "corrects" this input by 
incorporating the appropriate delay for each fraction of 
inputing links, h{r)d,T, with delay r. 

Expanding f(u>,6,t) in a Fourier series, we have 
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where = J Q 2 ^ f{uj 1 9 1 t)d9 is the time-independent 

oscillator frequency distribution. Following the method 
outlined in [13j |. we consider the dynamics of Eq.© on 
an invariant manifold in /-space: 



f n (LU,t) = [a(u,t)] n . 



(7) 



The macroscopic dynamics of a(u>,t) can be derived by 
substituting Eq.© into Eq.Q, leading to 



da/dt 



(A/2) 



= 0. 



(8) 



In the case when the oscillator frequency distribution 
g(u)) is Lorentzian, i.e., 



A/tt 



(w - uj ) 2 + A 2 



(9) 



and assuming suitable properties of the analytic continu- 
ation into complex uj of a(uj,t) (see Ref.[13j), Eq.© can 
be evaluated explicitly by contour integration with the 
contour closing at infinity in the lower half complex em- 
plane to give = J_ g(u>)a*(u),t)du) = a*(cJo~ iA,t). 
Thus Eq.© becomes 



r(t) 



a*(t — T)h(r)dT. 



(10) 



Furthermore, by setting lu = luq — iA in Eq.© we have 
j t a{t) + (A + two)o(t) + ~ (r(t)a(t) 2 - r* (t)) = 0, (11) 

where in both Eqs. lTOl) and (fTTj) the particular argument 
value to = LUQ — iA has been suppressed; i.e., a(wo — iA, t) 
is replaced by a(t). Equations (fT0|) and ([TT]) thus form 
a complete description for the dynamics on the invariant 
manifold when o(w) is Lorentzian. Recently a result 
has been obtained [14j that, when applied to our problem, 
establishes that all attractors of the full system, Eqs.© 
- ©, are also attractors of our reduced system, Egs-lfTO]) 
and (jTTj) , and vice versa. (The result of Ref. 1J] was 
previously strongly indicated by numerical experiments 
of Ref. 



Previous studies of the effect of delay on the Kuramoto 
system (Refs.[9l. [Tol 11 j) considered uniform delay on all 
the links, corresponding to h(r) = 5(r — T). Our goal 
is to uncover the effect of heterogeneity of delays along 
the network links. Accordingly, we consider that h{r) 
has some average value T with a spread about this value, 
and h(r) = for r < 0. A convenient class of functions 
for this purpose is 



Kt) = ^h n (1) , where ^(f) = A n f n e~^ f . ( 



12) 



Here, A n and (3 n are determined by the normalization 
conditions: J Q h n (f)df — 1 and J Q Th n (r)dT = 1, 
yielding 



A n = (n+l) n+1 /n\ , n =n + l. 



(13) 



For this family of distributions, we have that the standard 
deviation of r about its mean T is given by 

2U/2 



5t = (< T 2 > - < T > Z 



T/Vr 



1. 



(14) 



Thus, for n — > oo, we recover the case, h(r ) = 8{r — 
T), previously investigated in Refs.0, [l0, EH- And, by 
decreasing n, we can study the effect of increasing the 
relative spread St /T in the delay times. The dependence 
of h(r) on n is depicted in FigfT] 
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FIG. 1: Graphs of h(r) at T = 1 for n = 1, 20, 100. 

We can exploit the convolution form of ([5]) and turn 
it into a differential equation for r(t). Taking a Laplace 
Transform, we have, for the case of Lorentzian g(u>), 



f(s)=H( s )d*(s), 



(15) 



where r(s) and a(s) are the Laplace transform of r(t) and 
a(t) respectively, while 



F( S ) = [(T//3„) S + l]-(" +1 ), 



(16) 



is the Laplace transform of h(r). Our choice of the 
function class given by Eq. (TT2")) is motivated by the fact 
that it yields a particularly convenient Laplace transform 
and corresponding time-domain formulation. In partic- 
ular, transforming back to the time-domain by letting 
s -> d/dt, Ea.lfIB) yields 



[(T/(3 n )(d/dt) + l] n+1 r(t)=a*(t). 



(17) 
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Thus, we now have Eqs. (fTT|) and (fT7j) as our descrip- 
tion for the dynamics on the invariant manifold with 
heterogeneous link delays. Here, it is noteworthy that 
Eqs. lfTTj) and (fTT|) form a system of ordinary differential 
equations in comparison with the original system Eq. JT]) 
which comprises a very large number of time-delay dif- 
ferential equations. Note that for the case of uniform 
delay, h{r) — S(t — T), we take the limit n — ► oo, in 
which case Eq. (fT"5|) takes the form f(s) — e~ sT a*(s), 
yielding r(t) = a*(t — T), which, when substituted into 
Eq. tTTTj), g ives the time-delay differential equation for a(t) 
in Ref.Qj. 

A trivial exact solution to the system (fTTj) and (fT7]) 
is given by r(t) = a(t) = 0, which we refer to as the 
"incoherent state" [la]. Stability of the incoherent state 
can be studied by linearizing Eq. (fTT|) about the solution 
a(t) 



a e° 



and setting a(t) 

[kH(s)/2](s + iu! + A) 



from which we obtain 



1 



(18) 



The critical coupling k c at which a stable incoherent state 
solution becomes unstable as k increases through k c , cor- 
responds a solution to Eq. (fl8|) with Re(s) = 0. 

The solid curves in Fig[2] show results obtained from 
Eq. (|18|) with Lorentzian g(uj) for the critical coupling 
value k c versus T at different n's with parameters loq = 3 
and A = 1. For the case of uniform delays (n — ► oo), k c 
as a function of T exhibits the type of dependence found 
in Ref.[9j with characteristic "resonances". However, as 
the relative spread St/T is increased (n is decreased), we 
see that the resonant structure that applies for the case 
of zero spread is strongly modified. For example, even 
at the relatively small spread of St/T ps 0.1 (correspond- 
ing to n — 100), there is only one peak (at T ps 1) and 
one minimum (at T ps 2), with k c for T > 2 being very 
substantially higher than in the case of no spread. For 
St/T ps 0.302 (n — 10) the effect is even more severe, and 
the previous resonant structure is completely obliterated. 
For comparison, the dashed curves in Figl^show results 
for St/T ps 0.302 (upper) and (lower) when g{w) is 
Gaussian with the same peak value as for the Lorentzian 
distribution used to obtain the solid curves (l7j |. The 
Gaussian and Lorentzian results are similar, suggesting 
that the qualitative behavior does not depend strongly 
on details of g(ui). 

As reported in Ref.Q, bistable behavior can exist; i.e., 
a situation in which both incoherent and coherent states 
are stable. In Figs 3(a) and |3(b)| we show the hystere- 
sis loops obtained by numerical solution of Eqs. (|lip and 
(|17[) for n < oo and, for n = oo, where the n — oo result 
is obtained by solution of the delay equation obtained by 
inserting r(t) = a*(t— T) in (fTTj). Comparing Fig 3(a)| 
which is for T = 1, with Fig 3(b) which is for T = 3, 
we note the striking result that, for large T, hysteresis is 
sustained only with large enough spread in the delay dis- 
tribution, i.e., when n is small [e.g., for n = oo and T = 3 
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FIG. 2: Solid curves are plots of the critical value of 
k c versus T for Lorentzian g(uj) with &q ~ 3, A = 1, 
and n = 10, 100, 500, 1000, oo, corresponding to 8t/T ~ 
0.302, 0.1, 0.045, 0.032, 0. The dashed curves are for Gaussian 
g(ui) as described in the text. 



(Fig |3(b)| the bifurcation is supercritical and hysteresis 
is absent]. 
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FIG. 3: (Color online) Hysteresis loop at loq = 3, A = 1, (a) 
for T = 1, (b) for T = 3. 



Coherent oscillatory states can be obtained by substi- 
tuting the ansatz a(t) — a^e 1 , where ag and are real 
constants, into Eq. lTUj) and (flT)) . This gives 

+ lu ) + A] +{k/2) (alH*(iQ) - H(ift)) = 0, 

r(t) =a e- iUt H*(iQ). (19) 

As reported in both [9] and [l(|, for Ii(t) = 5(t - T), 
multiple branches of coherent state solutions are possible 
in EQ. (fT9"|) . Furthermore, we can employ Eqs. fTTj) and 
(fTTj) to study the stability of each coherent state by in- 
troducing a small perturbation 8a{t)e lQ,t to the coherent 
state solution in |[T9"]). with Sa(t) = K%e st + K 2 e s ' t . This 
yields the following equation for s: 



[s - ^H(s + iil) + A][s - |if(s - iSl) + A* 
= (kal/2) 2 H(s - iQ)H(s + id), 



(20) 



where A = IS.+i{ujQ-\-Vt)+ka^H (— iVt). Instability of each 
coherent state is then determined by whether there are 



solutions to (j2TJj) for s with positive real parts. In Fig 4(a) 
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we compare the theoretical results for \r\ calculated from 
Eqs. (fl9f and (|20| with simulation results based on Eq.|T|) 
with TV = 100 and 5t/T « 0.1 (n = 100) for the first two 
branches of coherent states with cj = 3, A = 0.1, T = 1. 
The solid (dashed) curves correspond to stable (unstable) 
coherent states. The Eq.((T]) simulation values reported 
in the figures represent time averages of these quantities 
computed after the solution has apparently settled into 
the coherent state. It is seen that there is good agree- 
ment between the theory and simulations using Eq.(fT]). 
In addition, on simulating these two branches of coher- 
ent states, we verified that the finding of Ref. [l(| that 
the basin of attraction is large for the first branch, but 
small for the second one, also holds with heterogeneous 
delays. 

Furthermore, the number of coherent attractors 
strongly depends on the spread in delay times. Figure 
|4(b)| shows the dependence of the number of coherent at- 
tractors on the relative delay spread St/T = (ri + l)~ 1//2 , 
with k = 40, wo = 0, T = 1, for two values of the fre- 
quency spread, A = 5 (dashed) and A = 10 (solid) (for 
which k c — 10 and 20, respectively). For both cases, 
it is seen that as the relative delay spread is increased 
((n + l) -1 / 2 is increased), the number of coherent at- 
tractors decreases. And there remains at least one such 
attractor when n approaches unity, while a parameter 
dependent maximum is attained when n — > oo, which we 
find is generally larger for smaller A and larger k [18j ]. 
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FIG. 4: (a) Magnitude of r for the first two branches of 
coherent states with parameter values: uo = 3, A = 0.1, 
n = 100, T — 1 for h(r). Solid lines give the theoretical values 
of the stable coherent states, dashed lines give the unstable 
coherent states, and asterisks give the simulation results, (b) 
Number of coherent attractors (number of solutions of Eq. (|19[l 
that are stable according to (|20|) ') versus St/T for the follow- 
ing parameters: k = 40, T = 1,cjo — 0; A = 10 for the solid 
line (k c = 20), A = 5 for the dashed line (k c = 10). 



In conclusion, in this paper we address, for the first 
time, the effect of heterogeneous delays on the dynamics 
of globally coupled phase oscillators. As compared to the 
case of uniform delay (Refs. @, M, [Tl| ) , we find that delay 
heterogeneity can have important consequences, among 
which are the following: (i) decrease in resonant structure 
of the dependence of k c on T (Fig|5|) ; (ii) increase of k c 
(Figl2); (iii) enhancement of hysteretic effects (Figs 3(a) 



and |3(b)[ ) ; (iv) reduction in the number of coherent at- 
tractors (Fig |4(bJ| ). Furthermore, we have introduced a 
framework for the study of delay heterogeneity that can 
be readily applied to a variety of extensions of the Ku- 
ramoto model, such as communities of oscillator popula- 
tions with different community dependent characteristics 
non-monotonic g(u>) [l5| , and periodic driving [2(| 
(see Ref.JlH for more examples). 
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